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Compact astrophysical objects that rotate rapidly may 
encounter the dynamical "bar instability." The bar-like de- 
formation induced by this rotational instability causes the 
object to become a potentially strong source of gravitational 
radiation. We have carried out a set of long-duration simula- 
tions of the bar instability with two Eulerian hydrodynamics 
codes. Our results indicate that the remnant of this insta- 
bility is a persistent bar-like structure that emits a long-lived 
gravitational radiation signal. 

PACS numbers: 04.30. Db, 04.40. Dg, 95.30.Lz, 97.60.-s 



I. INTRODUCTION 

The direct detection of gravitational radiation presents 
one of the greatest scientific challenges of our day. 
With interferometers such as LIGO, VIRGO, GEO, and 
TAMA [0 expected to be operating in the next few years, 
and a new generation of spherical resonant mass detectors 
under study |^,^ , the calculation of the signals expected 
from various astrophysical sources has a high priority. 
Accurate calculations of the waveforms are needed to en- 
able both the detection and identification of sources . 
In particular, short duration bursts are expected to be 
more difficult to detect than longer- lived signals. 

One interesting class of sources includes rapidly ro- 
tating compact objects that develop the rotationally- 
induced "bar instability" . This instability derives its 
name from the bar-like deformation it induces. The resul- 
tant object is potentially a strong source of gravitational 
radiation because of its highly nonaxisymmetric struc- 
ture. Examples of compact astrophysical objects that 
may rotate rapidly enough to encounter this instability 
include stellar cores that have expended their nuclear fuel 
and are prevented from undergoing further collapse by 
centrifugal forces [sj-p^ ; a neutron star spun up by accre- 
tion from a binary compani on ||ll|, [l2||; and the remnant 
of a compact binary merger |13| , ]l4[[ ! 

Such global rotational instabilities in fluids arise from 
nonaxisymmetric modes e^*""^, where m = 2 is known 
as the "bar mode" . It is convenient to parametrize them 
by 



/3 = r,c 



\w\ 



(1) 



where Tj-ot is the rotational kinetic energy and W is the 
gravitational potential energy In this paper, we 

focus on the dynamical bar instability, which is driven by 
Newtonian hydrodynamics and gravity, and is expected 
to be the fastest growing mode. It operates for fairly 
large values of the stability parameter P > f3d and de- 
velops on a timescale on the order of the rotation period 
of the object. For the uniform density, incompressible, 
uniformly rotating Maclaurin spheroids, (3d w 0.27. In 
the case of differentially rotating fluids with a polytropic 
equation of state, the m — 2 dynamical stability limit 
(3d ~ 0.27 has been numerically determined to be valid for 
initial angular momentum distributions that are similar 
to those of Maclaurin spheroids Jl6|-[T9t ; see also |^,^ . 
(We note that when (3 is greater than some critical value 
/3s < /3d, a secular instability can arise from dissipative 
processes such as gravitational radiation reaction and vis- 
cosity. When this instability arises, it develops on the 
timescale of the relevant dissipative mechanism, which 
can be several rotation periods or longer |jl^. In recent 
years, much work has also been carried out on various 
other modes in rotating relativistic stars as detectable 
sources of gravitational radiation; see |2^ for a review 
and references.) 

The first numerical simulations of the dynamical bar 
instability were carried out by Tohline, Durisen, and Mc- 
Collough (TDM; ||2^) in the context of star formation. 
Using a polytropic equation of state. 



Kp^ ^ Kp 
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with polytropic index N — 1.5, they evolved differentially 
rotating axisymmetric models with a 3-D Eulerian hydro- 
dynamics code, or hydrocode, in Newtonian gravity. In 
all models with initial /3 > 0.30, the m = 2 mode grew 
to nonlinear amplitudes and a two-armed spiral pattern 
was produced as mass and angular momentum were shed 
from the ends of the bar. Numerous other simulations 
have confirmed this basic scenario for the evolution of 
the bar mode into the nonlinear regime; see Sec. g for 
references and further discussion. 

More recently, these techniques have been extended 
to the context of rapidly rotating, compact objects in 
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the Newtonian regime, with the gravitational waves cal- 
culated in the quadrupole limit. This is a reasonable 
first approximation for an object, such as a centrifugally- 
hung stellar core with a density intermediate between 
that of white dwarfs and neutron stars, with initial 
mass M — 1.4M0 and radius R ^ 100 km, and hence 
GM/Rc^ <^ 0.02. Centrella and collaborators used both 
smooth particle hydrodynamics (SPH) and Eulerian fi- 
nite difference hydrodynamics to evolve the bar insta- 
bility in a model with N = 1.5 l24l|2l. In aU of their 



runs the gravitational wave signal was a relatively short 
duration burst lasting for several bar rotation periods, 
and the system evolved to a nearly axisymmetric central 
core surrounded by a flattened, disk-like halo. New p^ ] 
carried out a similar study, with an improved version of 
Tohline's Eulerian code ||2^. Her simulation employed a 
symmetry condition that only permitted the growth of 
even modes to; see Sec. This simulation produced a 
final state with a persistent bar-like core, which yielded 
a gravitational wave signal of much longer duration than 
that found by Centrella and collaborators. 

Given the requirements of reliable waveforms for the 
detection and identification of sources, it is important 
to resolve this issue of the late-time gravitational wave 
signal from the dynamical bar instability. To this end, 
we have carried out a set of long-duration runs using the 
two Eulerian codes employed by New and by Centrella 
in their earlier work, and have made a detailed study of 
the resulting models. In Section ^ we review previous 
numerical studies of the dynamical bar instability, high- 
lighting the various assumptions and restrictions used by 
different authors. The numerical techniques we used are 
discussed in Sec. |l|. In Sec. |^ we present our simu- 
lations and analyze the results. A discussion of these 
results follows in Sec. 0. Finally, the Appendix contains 
additional information about the two hydrocodes used in 
this work. 



II. PREVIOUS NUMERICAL STUDIES 

As mentioned above, the work of TDM |2^ set 
the stage for subsequent numerical calculations of the 
dynamical bar instability. Their initial models con- 
sisted of differentially rotating, axisymmetric equilibrium 
spheroids with a Maclaurin rotation law for the angular 
momentum distribution. The Maclaurin law produces 
rigid rotation when it is applied to an incompressible 
(N = 0) fluid; when it is used in a polytrope, it pro- 
duces differential rotation |l^. After small amplitude 
random perturbations were applied to the density, each 
model was evolved into the nonlinear regime using a 3-D 
Eulerian hydrocode with Newtonian gravity. This hy- 
drocode solved the mass continuity and Euler equations 
in cylindrical coordinates (n7, z,(^); the resulting evolu- 
tions were adiabatic and maintained the same polytropic 
equation of state, Eq. (H). 



TDM used equatorial symmetry and "7r-symmetry" in 
their simulations. Equatorial symmetry is a reflection 
symmetry through the equatorial plane z = 0. The tt- 
symmetry condition imposes periodic boundary condi- 
tions at angles ip = ir and ip — 27r; thus, physical vari- 
ables are the same in the interval < iy9 < tt as they 
are in n < ip < 27r. It is computationally advantageous 
to impose an equatorial- and/or tt— symmetry condition 
on such a simulation because, by doing so, half as many 
computational grid zones are required in order to achieve 
a given spatial resolution in the vertical and/or azimuthal 
coordinate directions, respectively. It was also physically 
reasonable for TDM to adopt both of these symmetry 
conditions because the eigenfunction (a pure, m — 2 bar- 
mode) to which their models were expected to be initially 
unstable had both equatorial and tt— symmetry p3| , p8[ . 
As we discuss more fully below, ultimately one would 
like to remove these computational constraints in order 
to test whether or not the physical outcome is sensitive 
to them. 

The first work to address the late-time development 
of the bar instability was published by Durisen, Gin- 
gold, Tohline, and Boss who ran simulations with 
(3 = 0.33 and /3 = 0.38 for N = 1.5 polytropes. They 
used three different 3-D hydrocodes: Tohline's Eulerian 
code as used in TDM, another Eulerian code developed 
in spherical coordinates by Boss, and an SPH code devel- 
oped by Gingold. Boss's code also enforced equatorial- 
and TT-symmetries but, being gridless, Gingold's SPH 
code imposed neither of these symmetries. However, the 
SPH simulations were limited to a very small number 
of particles, A'p = 2000. The results produced by these 
three separate simulation codes were qualitatively simi- 
lar. For example, at early times all simulations showed 
evidence of the development of a bar-like pattern in- 
stability, consistent with the results of TDM and the 
predictions of linear perturbative analysis |l5|j23| , |28| , pT| . 
Peturbative analysis says this instability is the result 
of the growth of a coherent bar-like wave that prop- 
agates around the system with a well-defined pattern 
speed, while material moves differentially through that 
pattern. At subsequent times in the simulations, the 
barmode distortion developed into a two-armed, trailing 
spiral pattern as described by TDM; when the spiral pat- 
tern reached a nonlinear amplitude, some relatively high 
specific angular momentum material was expelled in the 
equatorial plane of each system; and the primary struc- 
ture that remained at the end of each simulation was a 
dynamically stable, centrally condensed object exhibiting 
a value of /3 < /3d- But there were significant quantita- 
tive differences among the various evolutions presented 
by Durisen et al. For example, the simulations produced 
central remnants that had different total masses and ex- 
hibited different degrees of nonaxisymmetric distortion. 
This disagreement signified, in part, that the simulation 
techniques being used were rather primitive and, in part, 
that the available computing resources did not permit 
the simulations to be carried out with adequate spatial 
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resolution. 

Williams and Tohline subsequently carried out an in- 
vestigation of the dynamical barmode instability in mod- 
els with different polytropic indices. Using the TDM code 
with TT-symmetry and an improved azimuthal grid reso- 
lution, they first considered models with initial /? — 0.31 
and N = 0.8, 1.0, 1.3, 1.5, and 1.8, and focused their anal- 
ysis on the measurement of barmode growth rates and 
pat tern speeds in the linear-amplitude growth regime 
[|||. The runs with N = 0.8 and N = 1.8 were 
then extended to later times through the development 
of nonlinear-amplitude nonaxisymmetric structures and 
yielded a rotating triaxial central remnant |^ . Williams 
and Tohline noted that such a configuration would be of 
interest when viewed in the context of compact stellar 
objects because "its existence would presumably be dis- 
cernable from the spectrum of any emitted gravity wave 
radiation," but they did not derive such a spectrum from 
their models. 

Houser, Centrella, and Smith ||2^,^ were the next to 
carry out 3-D simulations of the dynamical bar instability 
for the case A'' = 1.5 and initial /3 « .30, this time in the 
context of rapidly rotating stellar cores. Using both an 
SPH and an Eulerian code, they considered the matter 
to be a perfect fluid with equation of state 

P=(r-l)pe, (3) 

where e is the specific internal energy, and solved an 
equation for the internal energy. Using artificial vis- 
cosity, they could account for the energy generation by 
shocks that occurs when the spiral arms form and de- 
flect the streamlines of the supersonically moving fluid. 
Routines were added to calculate the gravitational wave- 
forms and luminosities in the quadrupole approximation. 
The SPH code (developed from TREESPH; see |l|) im- 
posed no symmetry restrictions and was run with up to 
A^p = 32, 914 particles. Their Eulerian code, written in 
cylindrical coordinates, imposed symmetry through the 
equatorial plane but not 7r-symmetry |3^ . Overall, their 
simulations produced nearly axisymmetric central rem- 
nants at late times. 

Houser and Centrella carried out additional SPH 
simulations with N = 1.5, 1.0, and 0.5, and initial (3 ~ 
0.30 using improved initial models with A^p 16,000 
particles. As before, the N = 1.5 case resulted in an 
almost axisymmetric central remnant and a correspond- 
ingly short burst of gravitational radiation. The runs 
with N — 1.0 and N = 0.5 underwent additional episodes 
of spiral arm ejection, with the number of episodes in- 
creasing as decreased; such behavior was also observed 
by Williams and Tohline |3^]. This resulted in longer- 
lived nonaxisymmetric structure in the central remnants, 
accompanied by longer duration gravitational waveforms 
as the models grew stiffer (i.e., as A' was decreased). Note 
that the relatively small number of particles present in 
the SPH simulations of Centrella and collaborators, ac- 
companied by the velocity dispersion in their initial mod- 
els, may make it difficult for models with softer equations 



of state (larger A') to maintain long-lived nonaxisymmet- 
ric structures. 

New used an improved version of Tohline's code 
Hi to study the N = 1.5, (3 = 0.30 case. This code 
solves an energy equation and incorporates artificial vis- 
cosity to handle the shocks. She added a routine to 
calculate the gravitational radiation in the quadrupole 
limit. Her simulation, which imposed both equatorial 
and TT— symmetries, produced a persistent bar structure 
and a long-duration gravitational waveform. 

All of the studies mentioned above in this section 
started from initially axisymmetric models with the same 
radial distribution of specific angular momentum as in 
a Maclaurin spheroid. Pickett, Durisen, and Davis |2^] 
studied the instabilities that result in an A^ = 1.5 poly- 
trope, when the angular momentum distribution is var- 
ied. They used a (different) updated version of Tohline's 
code with neither equatorial plane symmetry nor tt- 
symmetry imposed; all their evolutions were adiabatic. 
Using the Maclaurin rotation law, they evolved a model 
with initial /? = 0.327 to late times, and obtained a bar- 
shaped central remnant. 

Recently Imamura, Durisen, and Pickett |Q have 
performed additional adiabatic simulations of dynami- 
cal instabilities in A^ = 1.5 and 2.5 polytropes with the 
Maclaurin rotation law, using the same hydrodynam- 
ics code used in ||2^. They focused on comparing the 
early phases of nonlinear mode growth in their runs with 
the predictions of quasi- linear approximations. Their 
high resolution simulations of A^ = 1.5 polytropes with 
/? = 0.304 and 0.327 both resulted in bar-like endstates. 

The properties and outcomes of the long duration bar 
mode runs with A^ = 1.5 and A'^ = 1.8 are summarized 
in Table || for convenience. All of the times reported in 
Table || are given in units of tc, where tc is defined as one 
central initial rotation period (cirp). When surveying 
the information catalogued in Table I, one should keep 
in mind that the identified "final" state has been reported 
at different evolutionary times in the various references. 
As this table emphasizes, over the fifteen years that have 
passed since the original Durisen et al. comparison pa- 
per ]29| , there remain significant quantitative differences 
among the results of various published simulations of the 
bar mode instability. In particular, as indicated by our 
comments under the "remarks" column, these previous 
simulations do not clearly indicate whether or not the 
end product of the evolution should be a central, steady- 
state structure that has a bar-like geometry. 

III. NUMERICAL TECHNIQUES 

A. Initial Axisymmetric Equilibria 

The new simulations of the dynamical bar instabil- 
ity presented here begin with rotating spheroidal models 
above the Maclaurin stability limit, (3 > /3d, constructed 
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in hydrostatic equilibrium. For fluids rotating about the 
2— axis with angular velocity fl = flizu), where zu is the 
distance from the rotation axis, the equations of motion 
reduce to the equation of hydrostatic equilibrium, 



-VP + V$ 
P 



hlW^ = 0, 



(4) 



where ^'(tn) = — J n'^{m)w dm is the centrifugal 
potential and ho is a constant. The gravitational poten- 
tial <i> is a solution to Poisson's equation, 



= AnGp. 



(5) 



The initial models for the runs discussed in this paper 
were constructed using Hachisu's Self-Consistent Field 
(HSCF; see also |2^) technique, which is a grid- 

based iterative method. To facilitate treatment of the 
boundary conditions, it uses the integral form of the hy- 
drostatic equilibrium condition, Eq. (j^. This gives 



H + <S> 



(6) 



where H = J p^^dP is the enthalpy of the fluid and C is 
a constant determined by the boundary conditions. The 
models are computed on a uniformly-spaced (ro, z) grid. 
The method requires an equation of state P = P{p)- For 
the polytropic relation in Eq. (||) , the enthalpy takes the 
form 



H = (l + N)Kp 



l/N 



(7) 



For purposes of comparison with earlier work, we fol- 
low Bodenheimer and Ostriker |^6| and adopt a specific 
angular momentum profile that is the same function of 
cylindrical mass as a Maclaurin spheroid, namely. 



= /lo [l - (1 - m{m)/M)^/^ 



(8) 



where M is the total mass of the system, m{u7) is the 
mass interior to cylindrical radius zu, the constant ho = 
5J/2M, and J is the total angular momentum. Hence, 
the centrifugal potential is, 



*(tn) 



1 - {l~m{m)/Mf/^ 



^dw. (9) 



Because the angular velocity is assumed initially to be 
only a function of vj, Lichtenstein's theorem implies that 
the configuration will have reflection symmetry through 
the equatorial plane . 

The HSCF method requires that two boundary points, 
A and B, on the surface of the model be selected . For 
spheroids, point A is set along vj at the equatorial radius, 
w{A) = tUE, and point B is set on the z— axis at the 
polar radius, z{B) ~ zp. The axis ratio zp/we is given 
as an input parameter; varying it produces equilibrium 
models with different values of (3. Points A and B set the 
boundary conditions for the solution of Eq. (^). Since p. 



P, and therefore H vanish on the surface of the polytropic 
fluid, we have 



H{A) = = C - - 



H{B) = = C - - hl-^{B). 



(10a) 



(10b) 



Once $ and ^' are known, Eqs. ( |10[) can be solved for the 
constants C and h^. 

The HSCF iteration process begins with an initial 
guess for p{-m,z), which also specifies the mass enclosed 
within each cylindrical radius m{m). Given p, the gravi- 
tational potential <i>(tn, z) is determined by solving Pois- 
son's equation, Eq. (||); see Ref. |3^] for details. Given 
771(07), the centrifugal potential ^'(ccj) is determined using 
Eq. (^. Then, C and /i§ are found from the boundary 
conditions, Eqs. (p^), and the enthalpy H is computed 
from Eq. (|^). Finally, a new density distribution is cal- 
culated from H by inverting Eq. (Q); this is used as input 
for the next iteration cycle. The process is repeated until 
fractional changes in C and and the maximum frac- 
tional change in H between two successive iteration steps 
are less than some threshold (in this work, 10^^). The 
virial error VE provides a measure of how well the en- 
ergy is balanced, and thus is indicative of the quality of 
the resulting equilibrium configuration. It is defined by 



VE = 2T + W + 3 / PdV, 



(11) 



where T is the total kinetic energy, and V is the volume 
of the model. The VEs for the models used here are 
- 10-3. 



B. 3-D Hydrodynamics Codes 

The simulations presented in this paper were carried 
out using two hydrocodes that employ Eulerian finite- 
differencing techniques to solve the equations of hydro- 
dynamics coupled to Newtonian gravity. The 2? (Drexel) 
hydrocode is the same one that was used by Smith, 
Houser, and Centrella [p5l in their studies of the bar in- 
stability, whereas the £ (LSU) hydrocode is the one that 
was used by New and Tohline In this section we 

briefly describe these codes, highlighting differences be- 
tween them that we believe to be most relevant to the 
analysis and discussion of our results. Further details on 
the V and C hydrocodes may be found in the Appendix. 

Both 3-D hydrocodes are written on uniform grids in 
cylindrical coordinates (nj, z,(/3). The 2? code assumes 
equatorial plane symmetry. The C hydrocode allows the 
use of both equatorial and 7r-symmetries, as discussed in 
Sec. Both codes handle the transport terms using sim- 
ilar monotonic advection schemes that are second-order 
accurate in space, and impose the same outflow boundary 
conditions on the edges of the grid. The V and C codes 
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both solve energy equations, using the perfect fluid rela- 
tion of Eq. (I) to calculate the gas pressure and artificial 
viscosity to handle shocks. Finally, both codes solve Pois- 
son's equation, Eq. (||), for the Newtonian gravitational 
potential $ with boundary conditions on the edges of the 
grid specified in terms of spherical harmonics. 

Eulerian codes typically require that the mass density 
in a grid zone never be zero, and thus fill the "vacuum" 
regions with a fluid having some small density, piow To 
facilitate the comparison of results in this paper, both 
codes impose essentially the same conditions in the "vac- 
cum" regions. The density is set to p = piow if the density 
drops below piow in a zone. The specific internal energy 
is similarly limited by e > eiow, where Eqs. (||) and (|^) 
give Clow = KPio^^/i^ ~ 1) and K is the polytropic con- 
stant of the initial model. In addition, the velocities in 
the low density zones must be limited to prevent them 
from becoming too large and thereby requiring very small 
timesteps through the Courant criterion . The veloc- 
ities are limited when p < pum = IQ^piow- Specifically, 
in cells where p < pum, and Vz are set to the value 
0.5cs,„iaa;, i/thcy exceed Cs^max- Here, Cs,max is the glob- 
ally maximum sound speed. Additionally, v^f, is set to 
zero in cells where p < pum and v^/zu > ftum- Here, 
fliim = where fig is the central rotation speed of 

the initial model. 

The codes do have a number of differences. The most 
important of these is that, as discussed in the Appendix, 
the hydrodynamical equations in the C code are writ- 
ten in flux-conservative form whereas in the T? code they 
are not. The accuracy of the C code is second-order 
in both space and time |l60|-|62|. However, while the T) 
code is spatially second-order accurate, the accuracy of 
its time evolution scheme is between first and second or- 
ders 1^,^. Finally, the V code was written in Fortran 
77 and optimized to run on Cray vector computers such 
as the C90 and T90 in single-processor mode. These 
machines use 64-bit words in single precision, which al- 
lows the use of very low densities in the vacuum region, 
typically piow = 10""'^°. The C code was developed for 
the parallel MasPar MP-1 computer, and was written in 
MasPar Fortran, which is MasPar's version of Fortran 90. 
The MasPar computers use 32-bit words in single preci- 
sion, and the £ code uses piow — 10^^ in the vacuum 
regions. 



a model with f3 = 0.298 > (3d- Models with two differ- 
ent resolutions were constructed for this work. The lower 
resolution model was computed on a grid with — 64 
radial and = 64 axial zones; its equatorial radius ex- 
tends out to zone j = 26 and its polar radius to zone 
k — 7. The higher resolution model was computed on 
a grid with N^^ = Nz = 128; its equatorial radius ex- 
tending out to zone j = 50 and its polar radius to zone 
fc = 12. (Note that = in radial zone j = 2 and z = Q 
in vertical zone k = 2.) The computations of both mod- 
els required 21 iterations. The 64 x 64 model had a virial 
error VE = 2.69 x lO'^, as defined in Eq. (|ll|); for the 
128 X 128 model, VE = 7.52 x lO"''. Density contours of 
this model in the x — z plane are shown in Fig. |l| and plots 
of the angular velocity ri(tn) (cf. Eq. (||)) and equatorial 
plane density profile p{zu, 0) are given in Fig. || We have 
normalized the equatorial radius and central density to 
unity. We = 1 and pc = 1. 

To prepare the initial data for evolution with a hy- 
drocode, the 2-D model is swept around in the azimuthal 
direction to create a 3-D axisymmetric model. Random 
perturbations are then imposed on the density to trigger 
the instability when the evolutions are run. Following 
TDM we set 



p{w, z. If) = peq[1 + ao/(ro, z, ip% 



(12) 



where peq is the density calculated with the HSCF code, 
flo = 10~^, and —1 < / < -|-1 is a random number. 

The perturbed initial models were then evolved with 
the hydrocode. Note that the models were initially cen- 
tered on the origin. All calculations used equatorial plane 
symmetry. Two simulations were performed with the 2? 
code. Model Dl had ^ Nz = = 64 zones while 
model D2 had twice the angular resolution, = 128; 
neither model used 7r-symmetry. Models Dl and D2 
both used the density piow = 10"^" in the "vacuum" 
regions. Three simulations were performed with the C 
code. Model LI was run without 7r-symmetry and used 
the same resolution as model D2; for comparison, model 
L2 was run with 7r-symmetry and = 64 and thus had 
the same angular zone size Aip as model LI. Finally, 
model L3 was run without 7r-symmetry and used twice 
the radial and axial resolution, — Nz — N^ — 128. 
The models run with the C code all used piow = 10~^. 
Some basic properties of these models are summarized in 
the first five columns of Table || for convenience. 



IV. ANALYSIS OF RESULTS 



A. Properties of the Models 

The initial axisymmetric equililbrium models used for 
the simulations presented in this paper were computed 
with the HSCF method [|5| as described in Sec. |III A . 
Specifying an axis ratio zpJwe = 0.208, polytropic index 
N — 1.5, and the Maclaurin rotation law, Eq. (||), yields 



B. Dynamical Evolutions 

All the models give similar results for the initial growth 
of the instability and its development into the nonlinear 
regime. An initially axisymmetric system develops bar- 
shaped structure as t he m = 2 mode grows to nonlinear 
amplitude (cf. Sec. IV C, [0). A pattern of trailing 



spiral arms is formed as mass and angular momentum 
are shed from the ends of the rotating bar. After the bar 
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reaches its maximum elongation, it recontracts somewhat 
towards a more axisymmetric shape. All the runs dis- 
played these basic characteristics, which are illustrated 
in Fig. 1^ using data from model L3. Time is measured in 
units of , where 

is the dynamical time for a sphere of mass M having the 
same radius as the initial equatorial radius of the model. 

Significant differences in the models arise in the next 
stage of the evolution, shown in Figs. ^ and |^. Since 
the models go unstable at slightly different times, they 
encounter the various phases of the instability at some- 
what different intervals; the frames in these figures are 
labeled with the time measured from the initial moment 
in each model. Frames (a-c) of Fig. ^ show that the cen- 
tral regions of model Dl appear to undergo a very slight 
re-expansion to a weak bar. After another ^1 — 2 bar ro- 
tation periods this weak bar disappears, leaving behind a 
nearly axisymmetric remnant that is somewhat displaced 
from the center of the grid. In model D2 the central bar 
re-expands more strongly, although not to the extent of 
its previous maximum elongation. After another ^1 — 2 
bar rotation periods, the model evolves toward a nearly 
axisymmetric remnant, which is offset from the center 
of the grid; see frames (d-f) of Fig. ^. Similar behavior 
is seen in model LI, although the second bar elongation 
phase is stronger and lasts somewhat longer; representa- 
tive density contours are shown in frames (a-c) of Fig. ||. 
The higher resolution model L3, shown in frames (d-f) 
of Fig. 0, also re-expanded to a fairly strong bar, which 
underwent two additional episodes of contraction and ex- 
pansion. The remnant in model L3 remained bar-like in 
shape for > 4 bar rotation periods, but eventually it also 
decayed and settled into a nearly axisymmetric remnant 
as it moved off the center of the grid. In contrast, model 
L2 (which was run with 7r-symmetry) retained a fairly 
strong bar for many (> 8) bar rotation periods. This 
model was run for a longer time than any of the others 
and experienced ~ 5 episodes of contraction and expan- 
sion. At the end of the simulation, model L2 still had 
a strong bar centered on the origin; see frames (g-i) of 
Fig. |. 

As is illustrated in Fig. ^, for all five model evolutions 
each episode of expansion and contraction of the bar is 
mirrored in the time-evolution of the stability parameter 
/9 — Trot/\W\. As the bar develops and expands, /3 drops 
and reaches a local minimum when the bar is at its max- 
imum elongation. Then /3 rises to a local maximum as 
the central regions recontract. This behavior occurs be- 
cause Trot oc lu!^ oc J^//. Hence, as the bar elongates its 
moment of inertia / increases which, assuming its angu- 
lar momentum is approximately conserved, reduces the 
rotational kinetic energy Trot. Each subsequent episode 
of bar re-expansion can also be associated with a local 



minimum in the global parameter f3. 

As has been illustrated by Figs. ^ and ^, in each sim- 
ulation except model L2, the "final" remnant appears to 
have moved off of the center of the grid. In each case 
this displacement is associated with measurable motion 
of the center-of-mass of the system (such center-of-mass 
motion was also seen in some of the simulations in ) . 
Fig. shows the position of the system center-of-mass 
^CM = [-'^CM + ^cm] ^^'^ ^ ^ function of time for each of 
the four affected models. (Note that Zcm = at all times 
in all of our simulations since they all employ equatorial 
plane symmetry.) In runs Dl and D2, the center-of-mass 
begins to move rather abruptly at a time t/tD ~ 20, 
then after moving to a location log(i?cM) ^ —0.5 (Dl) 
and ~ —0.7 (D2), the center-of-mass motion slows con- 
siderably. Because both models have radial zones of size 
Aw 0.04, this location corresponds to approximately 
8 and 5 radial grid zones, respectively, or in terms of 
the initial equatorial radius of the model, ~ j'^e and 
^ ^taJe, respectively. For models LI and L3, the center- 
of-mass motion does not appear to level off, as seen in 
Fig. 0(b) . Notice that the center-of-mass moves farthest 
from the center of the grid in model LI, which has the 
same resolution as D2, and that the onset of this mo- 
tion is significantly delayed when the radial resolution is 
doubled in model L3. 

Since these simulations all have Newtonian gravita- 
tional fields and the systems are assumed to be isolated 
from the external environment, in each case the center-of- 
mass should remain fixed at the origin if the total mass of 
the system remains confined within the boundaries of the 
computational grid {i.e., if the system remains isolated 
from its environment) and if the equations governing the 
dynamics of an isolated system are being properly in- 
tegrated forward in time. Although some (< 5%) mass 
(and associated linear and angular momentum) does flow 
radially off the grid, this mass loss does not appear to be 
large or asymmetric enough to account for the center-of- 
mass motion. We have verified this by running a simu- 
lation (not presented in detail in this manuscript) with 
the same resolution as model LI but with an expanded 
grid {N^ ~ Nz ~ 128). The center of mass motion of the 
model evolved on the expanded grid was not significantly 
different from that present in the run on the smaller grid 
[N^ = Nz = 64); the mass lost from the expanded grid 
was less than 1%. 

Instead, we suspect that the center of mass motion 
arises from numerical errors. In "flux-conservative" 
finite-differencing schemes, such as the one employed here 
in the C code (see the Appendix), the advection term is 
handled in such a way that, for example, momentum and 
mass are guaranteed to be globally conserved if the dy- 
namical equations contain no source or sink terms. How- 
ever, such algorithms are not explicitly designed to pre- 
serve the position of the center-of-mass of the system 
and source terms due to gradients in the pressure and 
gravitational potential do naturally arise (see, for exam- 
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pie, the right-hand-sides of equations (B1)-(B3)). So dis- 
crepancies that inevitably will arise between the finite- 
difference representation of the dynamical equations and 
their exact differential counterparts can lead to center- 
of-mass motion that is unphysical. (An exception to 
this is the case where 7r-symmetry is explicitly enforced. 
With TT-symmetry imposed, odd azimuthal modes can- 
not grow and, in particular, no center-of-niass motion 
can develop.) Once such motion begins in either the V 
or C code, it evidently has a tendency to amplify rather 
than to damp. 

It is also instructive to examine the conservation of 
total angular momentum J. All models show some de- 
gree of angular momentum loss as the bar expands into 
the "vacuum" regions of the grid. To quantify this, we 
consider the loss of J in each model at the time that 
(3 reaches its first local minimum, and thus at the time 
the bar reaches its point of maximum expansion; see, for 
example. Fig. ^. (Shortly after this time, the models typ- 
ically lose mass as the spiral arms expand and material 
flows off the grid.) In general. Table || shows that models 
run with the V code lose more angular momentum dur- 
ing this stage than those run with the C code. Additonal 
tests with the V code showed that this loss of angular 
momentum increases as piow increases p5| . Overall, we 
attribute the better angular momentum conservation ob- 
tained with the £ code to the fact that it is written in 
an explicitly flux-conservative form whereas the V code 
is not; see the Appendix. 



C. Analysis of Fourier Components 

We can quantify the development of the dynamical in- 
stability shown in the previous section by examining var- 
ious Fourier components in the density distribution. To 
this end, the density in a ring of fixed vu and z can be 
written using the complex azimuthal Fourier decomposi- 
tion 



+00 



(14) 



where the amplitudes Cm of the various components m 
are defined by ppCl] 



Cm(tn, z) 



1 

2^ 



pi'!j7,z,ip)e-'""^d(p. (15) 



We shall also find it useful to define the normalized am- 
plitude 



l^ml = |an|/|Co| 



(16) 



where Co(w, z) = pi'cu, z) is the mean density in the ring. 
The phase angle of the m"^ component is defined by 



When nonaxisymmetric structure propagating in the f- 
direction develops into a global mode we can write the 
phase angle as 

4>m = CTrnt, (18) 

where cr,„ is called the eigenfrequency of the m"^ mode. 
The pattern speed is then 



HVi(n7, z) = — — = — 
m at m 



(19) 



(?!)m(ti7, z) = tan ^ [Im(A„)/Re(An)] 



(17) 



and the pattern period is — 2tt/W„i ]2^ , p0| . For the 
m = 2 mode, <J2 is twice the angular velocity of the bar 
and, hence, the rotation period of the bar is T2 = 47r/cr2. 

In practice, we implement Eq. ( [l5| ) in the codes by 
summing up the contributions over the azimuthal coor- 
dinate (0 < iy9 < 27r) in rings of width Atn centered on 
the origin at various values of vj in the equatorial plane 
z = 0. By examining the amplitudes of the Fourier com- 
ponents \Am\ at various distances from the rotation axis, 
we can determine when global modes arise in the mod- 
els. Since the rings are always centered on the origin, 
care must be taken when interpreting the results once 
the center-of-mass motion becomes significant. 

Figure || shows the growth of the amplitudes |A,„| for 
the first four Fourier components, m = 1,2,3,4. These 
amplitudes were calculated in the equatorial plane in a 
ring of radius -cu = 0.354, which corresponds to radial 
zone j — 10, for models Dl, D2, LI, and L2. For model 
L3, the amplitudes were calculated in a ring of radius 
m = 0.344, which is radial zone j — 18. Similar plots 
were obtained for rings at other values of nj within the 
central regions, indicating that these are all global modes. 

Notice that the exponential growth of the m = 2 bar 
mode (thick solid line) dominates all evolutions, as ex- 
pected from visual inspection of the density contours 
shown in Figs. ^ - [H In addition, all models show a clear 
TO = 4 mode (thin solid line) that begins growing expo- 
nentially once the bar mode is well into its exponential 
growth regime. The to = 2 and m = 4 modes both reach 
their peak amplitudes at about the same time in all runs. 
The odd modes appear at later times in all models except 
L2, in which 7r-symmetry was imposed (which prevents 
the development of odd modes; see Sec. O). The to = 1, 
or translational, mode (dotted line) begins growing after 
the TO = 4 mode. Comparing Figs. ^ and || shows that 
the TO = 1 mode grows as the center of mass motion in- 
creases, as expected. The to = 3, or pear, mode is shown 
by the dashed line. In models LI and L3, the m = 3 
mode is the last one to grow. In models Dl and D2, the 
TO = 3 mode grows at early times, but then lags behind 
the m = 1 mode. 

As mentioned in Sec. ||, analysis of Maclaurin sphe- 
riods suggests that models near the dynamical stability 
limit are only unstable to the m = 2 mode (higher order, 
even harmonics of the to = 2 mode may arise in the sub- 
sequent evolution). Models near this stability limit are 
not physically susceptible to the growth of odd modes. 
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This has been confirmed again recently by the perturba- 
tive analysis of Toman et al. (TIPD; In fact, TIPD 

demonstrate that TV = 1.5 polytropes with the Maclau- 
rin rotation law [Eq. ^ are only unstable to the m — 3 
mode when /3 is ^ 0.32. In general all our models (for 
which /3 — 0.298, initially) conform to this expectation, 
with the odd modes (which here are numerical artifacts) 
growing earlier in the models with lower resolution. In 
the L3 simulation, which has the highest resolution, the 
m = 1 and m = 3 modes develop at late times and are 
cleanly separated from the m = 2 and m = 4 modes. 

The growth rate d\n\Am\/dt for the to = 2 and 
TO = 4 modes can be determined by fitting a straight 
line through the data points in Fig. || during the time in- 
tervals in which In \Am \ is growing linearly with time. We 
find that all models have approximately the same growth 
rate for these modes, as seen in Table 01. To find am for 
these modes we plot cos 0m versus time and use a trigono- 
metric fitting routine to get cf. [p3|| . The function 
cos 4>m was used because (pm itself is multi- valued due to 
the tan~^ in Eq. ([1^). The fit was performed over the 
same time interval, chosen "by eye" , used to calculate the 
corresponding mode growth rate. Table III shows that 



models Dl and D2 have smaller eigenfrequencies and pat- 
tern speeds than those run on the C code. However, in 
all models we find that the pattern speeds for the to = 2 
and TO = 4 modes are nearly the same, 14^2 ^ W4. This 
confirms that the to = 4 mode is a harmonic of the bar 
mode and not an independent mode, as mentioned above 
and first pointed out by Williams and Tohline 

It is interesting to compare the data from our runs 
with the results of TIPD, who used a perturbative, lin- 
earized Eulerian scheme to calculate mode growth rates 
and eigenfrequencies of differentially rotating polytropes. 
As discussed by TIPD, this method produces much more 
precise results than traditional Lagrangian normal mode 
analysis (including the tensor virial approach [^ , 
which has been demonstrated to be inappropriate for dif- 
ferentially rotating polytropes) . TIPD calculated the bar 
mode growth rate and eigenfrequency for an axisymmet- 
ric N = 1.5 polytrope with the Maclaurin rotation law, 
Eq. (P), and f3 — 0.300 using their perturbative Eule- 



rian method. They found dhn.\A2\/dt 
,-1 



(J2 



1.99 to', where we have converted from their units. 



TIPD state that the uncertainties in their results are on 
the order of a few percent (excluding any systematic er- 
rors that may be present). Comparison with the data 
in Table [II shows that the numerical models all have 



growth rates in good agreement with that of TIPD. The 
eigenfrequencies a2 for the models run with the C code 
are also in very good agreement with the TIPD results, 
while those for Dl and D2 are ^ 25% smaller. The rea- 
son for the discrepancies in the eigenfrequencies of the T) 
code runs is unknown. 



D. Gravitational Radiation 

We use the quadrupole approximation, which is valid 
for nearly Newtonian sources |^^, to calculate the grav- 
itational radiation produced in our models. Since the 
gravitational field in both codes is strictly Newtonian, 
we compute only the gravitational radiation produced 
and do not include the effects of radiation reaction. The 
reduced or traceless quadrupole moment of the source is 



d'r, 



(20) 



where i, j = 1,2,3 are spatial indices and r — {x^ + + 
z^Y^"^ is the distance to the source. For an observer lo- 
cated on the axis at0 = O,(/3 = Oofa spherical coordinate 
system with its origin located at the center of mass of the 
source, the amplitude of the gravitational waves for the 
two polarization states becomes simply [p5p^] 



r 

h - ^ 2 • • 



(21) 
(22) 



where an overdot indicates a time derivative d/dt. 

A straightforward application of Eqs. (|^ - ( |2l| ) in an 
Eulerian code involves calculating /.y directly by sum- 
ming over the grid, and then taking the time derivatives 
numerically. However, such successive application of nu- 
merical time derivatives generally introduces spurious 
noise into the waveforms, especially when the timestep 
is allowed to change from cycle to cycle. To reduce this 
problem, we have used the partially integrated versions 
of the standard quadrupole formula developed by Finn 
and Evans The V code uses the first moment of 

momentum formula, QFl, which allows the calculation 
of iij directly from quantities available in the code [ p2[ . 
This gives 



f — 2 



P [v 



d^T 



0.532 tj-, and where 



(23) 



(24) 



Of course, another time derivative is still required to ob- 
tain Iij . When this is taken numerically, the resulting 
waveform amplitudes h-^ and hx can still be dominated 
by noise. To cure this problem, we pass Iij through a 
filter to smooth it before calculating the waveforms; see 
Ref. H for details. 

The C code does not use numerical time derivatives to 
compute Iij. Instead, the equations of motion are used 
in conjuction with Eq. (|2^ ) to compute f y directly from 
quantities available in the code. This gives 
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ilm = J p[2viVm - -^Rv'^ViSlm ~ Xm^ 

--a;*Vi$(5,™ + (A.V. terms)]d^x. (25) 
3 

Here, the summation convention is implied on repeated 
up and down indices and the "A. V. terms" contain con- 
tributions to the stress tensor from the artificial viscos- 
ity; see Refs. and for details. The expression 
for Ilm in Eq. (E5h yields smooth waveforms which do 
not require filtering. Note that both the C and V codes 
compute lij with respect to the origin of the coordinate 
system. Hence when the center-of-mass motion becomes 
significant, the waveforms computed are no longer pre- 
cisely correct. However, the center-of-mass motion itself 
is a numerical artifact; thus, upon its development the 
simulation as a whole becomes physically unreliable. 

Fi gur e ^ shows the gravitational waveform hj^ given in 
Eq. as a function of time for all models. At early 
times, the waveforms are all similar, as the initial ex- 
pansion of the bar gives rise to gravitational radiation. 
Comparing Fig. ^ with Fig. we see that the maximum 
amplitude of occurs at the same time as the first 
local minimim in /3, which marks the maximum expan- 
sion of the bar. The amplitude of /i+ then dips as the 
central bar recontracts and /3 rises again to a local max- 
imum. When the bar re-expands, the amplitude of /i+ 
rises again; cf. Figs. ^ and ||. In models Dl, D2, and LI, 
the central remnant soon becomes nearly axisymmetric 
and /i+ decays rapidly. In L3, the bar persists for > 4 bar 
rotation periods and undergoes two additional episodes 
of expansion and contraction, producing a longer-lived 
gravitational wave signal. Model L2 maintains a strong 
bar with multiple expansions and contractions, and thus 
a strong gravitational wave signal, throughout its evolu- 
tion. 

In Fig. |l0| we plot the quantity 

Ko.m^{h\ + hlY/^. (26) 

Note that the absolute scale of the t axis in this figure 
is not labeled because the curves for the three runs have 
been shifted horizontally in order to line up their initial 
peaks. All models show a strong initial peak in /inorm that 
coincides with the maximum expansion of the bar. The 
secondary peaks in /inorm correspond to the secondary 
expansions of the bar and the additional local minima in 
(i shown in Fig. ^. A plot of /inorm is instructive because 
its time variation is not complicated by the periodic ro- 
tation of the bar. Thus ft.norm(i) refiects how the "mean" 
properties of the system (e.g., the moment of inertia) 
change with time. The /inorm curve should be perfectly 
fiat once (and if) the remnant settles down into a steady- 
state structure; any slight downward slope will provide a 
measure of long-term secular changes. 

Ultimately, every hydrodynamics code should pro- 
duce the same, correct, plot of /inorm(i) for this simu- 
lation. That is, the amplitude and frequencies present in 
^norm(0 should be quantitatively reproduceable. Fig. |l^ 



thus shows to what degree our (three best) simulations 
are converging toward the same answer and thereby be- 
gins to establish the nature of the "correct" result. 

V. DISCUSSION AND CONCLUSIONS 

We have carried out simulations of dynamical instabil- 
ity in a rapidly rotating = 1.5 polytrope using two dif- 
ferent Eulerian hydrocodes and several different resolu- 
tions. The rapidly rotating polytropic initial models used 
were constructed with the Maclaurin rotation law and 
had a ratio of kinetic to gravitational energy (3 ^ 0.3. All 
models evolved by both codes agree on the following basic 
properties of the early nonlinear development of the in- 
stability, (i) The m — 2 mode dominates the evolutions, 
producing a central rotating bar which sheds mass and 
angular momentum at its ends to produce a spiral arm 
pattern. Once the bar reaches its point of maximum elon- 
gation, it contracts and then re-expands, (ii) The m = 4 
mode is the next one to reach nonlinear amplitudes. (Hi) 
The growth rates for the m = 2 and m = 4 modes are 
d\n\A2\/dt « 0.55^0^ and d\n\A4\/dt w l.Ot^^, respec- 
tively. The pattern speeds W2 ~ W4, indicating that the 
TO = 4 mode is a harmonic of the bar mode, (iv) The in- 
stability produces a gravitational wave signal with maxi- 
mum amplitude [Re{c^ /GM)'^] rhK, 0.6 for an observer 
on the axis dX 6 — ip = i) oi a. spherical coordinate 
system centered on the source. 

The models also exhibit some differences. In particu- 
lar, simulations run with the 2? code have smaller values 
of the eigenfrequencies a2 and (74 , show weaker secondary 
bars, and lose more angular momentum during the ini- 
tial bar expansion than those run with the C code. It 
appears that most of these differences can be attributed 
to the lower order time differencing and the lack of con- 
sistent flux-conservative differencing in the V code. 

Overall the simulations presented here, and those car- 
ried out by previous authors, agree on the qualitative 
nature and many quantitative aspects of the initial de- 
velopment of the bar structure. However, as detailed 
in §1 and §11, such agreement has not been universally 
present among simulations that follow the long duration 
evolution of this instability. Such lengthy evolutions are 
nontrivial for hydrodynamics codes as they may allow 
any numerical inaccuracies present to grow to the point 
where they signficantly infiuence the simulations. As de- 
scribed below, the late growth of odd modes in our bar 
mode evolutions is an example of the numerical difficul- 
ties that can arise in extended simulations. 

Linear analysis indicates that the N = 1.5, (3 = 0.3 
models used in our simulations are initially unstable to 
the m — 2 bar mode only, and not to any odd modes 
pl| (even harmonics of the m=2 mode may subsequently 
develop). Thus the odd modes that develop in all but 
one (L2) of the simulations presented here are numerical 
artifacts arising from shortcomings in the finite-difference 
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techniques utiHzed in the T) and C hydrocodes (see Sec. 
IV.B and also 1^). 

Once these artificial odd modes reach nonlinear am- 
plitudes, the physical reliability of the simulations is de- 
graded. In particular, the growth of the m = 1 mode 
is tied to the development of center-of-mass motion in 
the simulations we ran without 7r-symmetry. As the 
center-of-mass moves away from the origin of the cylindri- 
cal grid, the finite-differencing of the curvilinear form of 
the hydrodynamics equations becomes asymmetric. This 
causes the accuracy of the evolution to deteriorate. The 
growth of the center-of-mass motion appears directly re- 
lated to the decay of the bar-like structure of the system; 
the bar decays when the center-of-mass motion becomes 
significant. 

Because the center-of-mass motion is unphysical, we 
believe the decay of the bar structure is unphysical as 
well. This conclusion is substantiated by run L3, which 
had twice the radial (and axial) resolution of run LI. The 
onset of the spurious center of mass motion was signifi- 
cantly delayed in L3 and the central remnant maintained 
its bar-like structure for a correspondingly longer time. 
The growth of odd modes was also delayed. Thus as the 
resolution is increased, the L code evolutions converge to- 
wards the predictions of linear analysis. (Unfortunately, 
we could not repeat this experiment with the T) code due 
to a lack of computational resources.) 

Thus it is our belief that a simulation of this instability 
that did not develop the nonphysical center of mass mo- 
tion (e.g., one performed with very fine radial resolution), 
would produce a long-lived nonaxisymmetric structure. 
Recall that this is the result of the model L2 simulation, 
which was run with 7r-symmetry. That symmetry condi- 
tion prevents the growth of odd modes and thus prohibits 
the development of center of mass motion. Hence, over- 
all L2 is the most physically relevant of the simulations 
presented here. 

We believe our results demonstrate that the physically 
acurate outcome of the rotational instability in the ob- 
ject studied here, is a persistent bar with an accompa- 
nying long-lived gravitational wave signal. This dynami- 
cally stable configuration can be viewed as a compressible 
analog of a Riemann ellipsoid; efforts to understand the 
detailed structural properties of such configurations are 
underway ]4^ . The nonaxisymmetric structure of the 
remnant will decay on a secular timescale due to viscous 
dissipation and/or gravitational radiation emission. The 
gravitational radiation timescale is likely to be shorter 
than the viscous timescale for sufficiently compact ob- 
jects 1 47 4^. Note that (3 will also decrease as a result of 
this secular evolution. The system will continue to evolve 
until it reaches a configuration that is secularily stable. 

A number of factors including the presence of an enve- 
lope surrounding the rotating object, the variation of the 
rotation law and the equation of state, and the influence 
of general relativity could potentially affect the outcome 
of the instability. Such effects should be the subject of 
future study. 
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APPENDIX A: 

The V and C codes solve the equations of hydrody- 
namics, which govern the structure and evolution of a 
fluid, in cylindrical coordinates. These equations include 
the continuity equation. 



dp 

at 



+ V • (pv) = 0; 



the three components of Euler's equation. 



dS 
dt 



V • {Sv) = - 



dP 
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dw dvj 
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— + V-{Tv)^-—^p—, 
at oz oz 
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dP 

dip 



5$ 
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(Al) 



(A2) 



(A3) 



(A4) 



and Poisson's equation, Eq. (||). Here, v is the fluid ve- 
locity with components (wro, I'z, tiip) in the [w,z,ip) di- 
rections. The quantities S = pv^, T = pv^, A = pwv^ 
are the radial, vertical, and angular momentum densities, 
respectively. 

The codes solve slightly different forms of the energy 
equation. The V code evolves the speciflc internal energy 
e: 
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dt 



1 d{wpev^) d{pevz) 1 d{pev^) 



-P 



dzu 
1 dijuv^ 
70 dvj 



dz 
) dv 

' A ■ 

dz 



zu dip 
1 dv^ 
zu dip 



The C code evolves the internal energy density e: 



de^ 
dt 



V • (e 



0; 



(A6) 



In both codes, the pressure P is obtained from the perfect 
fluid equation of state, Eq. (^. 

Both the V and C codes use artificial viscosity to 
smooth out sharp discontinuities that may arise if shocks 
are present in a simulation. See |^,^,^ for details. 

The following subsections contain further details about 
the V and C hydrocodes. 



1. The V Hydrocode 

The V hydrocode was developed by Clancy, Smith, 
and Centrella |^,|5^,^. It is written in cylindrical co- 
ordinates {zu, z, ip) with reflection symmetry through the 
equatorial plane z = 0. The original version allowed 
nonuniform radial and axial grids and was used to carry 
out the Eulerian runs in Ref. ; for the simulations de- 
scribed in this paper, the code was restricted to uniform 
grids. This code was written in Fortran 77 and optimized 
for Cray vector computers; it currently runs on the Cray 
T90. 

The actual form of the hydrodynamics equations ( |A1[ )- 

, with the 



(A4) used in the V code is given in Ref. 
exception that Eq. (A4) takes the form 



djpJ) 
dt 



1 d{pjv^w) d{pjv^) 1 d{pjv^) 



dvj 



dz 



UJ 

dP 

dip 



dip 
^ dp ' 



(A7) 



where J — zuv^ is the specific angular momentum. Note 
that the equations the T> code solves are not written 
in flux-conservative form. In the discrete form of these 
equations, the scalar quantities p, e, <I>, and P are de- 
fined at cell centers and at integral timesteps. The T> 
code actually evolves the velocities, which are defined 
on the faces between cells and at half-integral timesteps, 
located halfway between the integral timesteps. The ve- 
locities are face-centered in the coordinate along which 
they are directed; for example, Vz is defined at the center 
of the grid zone faces normal to the z— axis [ p2[ . 

The T) code uses operator splitting to evolve the 
discre te versions of the hydrodynamical equations, 
Eqs. @ - ®, forward in time |3|,|5|. The accu- 
racy of this time integration method is between first and 
second orders. 

The source step is carried out first. This begins by 
holding p constant and updating the velocities due to 



the pressure gradient, gr avit ati onal force, and centrifu- 
gal force terms in Eqs. (A2) - (A4) using centered dif- 
ferences; note that in the source step we advance the 
azimuthal velocity component instead of the specific 
angular momentum J in Eq. (A7). Using these updated 
values, the artificial viscosity terms are applied to ad- 
vance the velocities and e. These new values are then 
used to update the energy due to the compressional or 
"PdV" terms. 

We next carry out the transport step to evolve p, pe, 
Wro, Vzt and J due to the advection of fiuid from one 
cell to the next. The transport step consists of three 
advection sweeps, one in each of the three coordinate 
directions. We use a monotonic advection scheme de- 
veloped by LeBlanc |^,^ that is second-order accurate 
in space to calculate the fluxes in each direction. On 
each cycle, we vary the order in which the advection 
sweeps are carried out to avoid setting up a preference 
for any one direction; the order changes on each suc- 
cessive cycle as all six permutations are exhausted. On 
each sweep, the same mass flux used to advect the den- 
sity i n Eq . ( [Al|) is employed to advect u^, w^, and J in 
Eqs. ( [A^ ) - ( [A4| ) |]5^-|55t. During the transport step, the 
density is held c ons tant; thus, d^pv^^jdt i s wr itt en a s 
pdv^/dt in Eq. (^), and similarly for Eqs. Q - (^). 
After updating the advection terms on each cycle, a mo- 
mentum conservation is applied with the new density to 
update the velocities. The equation of state, Eq. (g), is 
then used to calculate a new value of the pressure. 

Once the hydrodynamical equations have been ad- 
vanced, the Newtonian gravitational potential $ is cal- 
culated by solving Poisson's equation, Eq. (||), using the 
updated density. The boundary conditions at the edge of 
the grid are specifled using a spherical multipole expan- 
sion. The discretization yields a large, sparse, banded 
matrix equation which we solve using a preconditioned 
conjugate gradient method with diagonal scaling |^,^. 



2. The C Hydrocode 

The £ hydrocode was originally developed by Tohline 
1^ , ^ , and has been reflned and updated with collab- 
orators and students. The modern version of the code 
is fully second order accurate in both space and time 
p7| , p9[ . The parallel version of the code that we use here 
was originally developed for the MasPar MP-1 computer 
and was written in MasPar Fortran, which is MasPar's 
version of Fortran 90; see The C code uses uni- 

formly spaced grids in cylindrical coordinates {zu,z,p). 
The code allows the use of reflection symmetry through 
the equatorial plane z = and tt— symmetry in the az- 
imuthal direction; cf. Sec. Ill 



The fluid equations, Eqs. (p\l 



|) and ( [Aq ) , are writ- 
ten in flux-conservative form |55t]. When they are dis- 
cretized on the uniform cylindrical grid, the density p, the 
angular momentum density A, and the gravitational po- 



ll 



tential $ are defined at cell centers. The radial and ver- 
tical velocities (v^,Vz) and momentum densities (iS,T) 
are defined at cell vertices or nodes. The source terms 
on the right-hand sides of Eqs. (A2) - (A4) are approxi- 
mated using standard second-order centered differences. 
The flux or divergence terms are written as a summation 
over the six faces of a cylindrical grid zone , 



1 ^ 



(A8) 



Here V is the volume of the cylindrical grid cell, Ai is the 
area of a particular face, and {Xv)i is the product of the 
quantity X e {p,S,T,A,e^^^) and the corresponding 
velocity component at the face i (i.e., v is the velocity 
normal to the i*'' face). These terms are updated using 
a monotonic interpolation scheme developed by van Leer 
|l60| that is second-order accurate in space. 

When the system is evolved forward in time, the phys- 
ical variables X S (p, S, T, A, e^^^) are updated by ap- 
plying the source terms and the flux terms in different 
steps. Second-order accuracy in time is obtained via a 
Lax-Wendroff scheme that uses velocity values in the flux 
terms (A_8) that are centered at time t + At/ 2 |^,^. 
To accomplish this, the source terms are applied for a 
half timestep At/2 and the updated values X' are saved. 
The flux terms are then applied for At/2 with the up- 
dated values X to obtain velocities at time t + At/ 2. 
With these new velocities, fluxing is performed for a full 
timestep on the saved quantities X' . An additional half 
timestep of sourcing is then performed. 

Poisson's equation, Eq. (||), is solved for the gravita- 
tional potential $ using the ADI method |^^, after the 
fluxes have been applied for the whole timestep At. 
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Fig.1 
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to 

FIG. 1. Density contours of the initial model with resolu- 
tion Nt^ — Nz ~ 128 are shown in the x — z plane. The 
maximum density is located at the center and is normalized 
to unity. The density contours are at levels of 0.5, 0.05, 0.005, 
and 0.0005. 
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Fig. 2 




03/03 E 

FIG. 2. The normalized angular velocity fl{zu)/Qc and 
equatorial plane density p{z = 0)/pc in the initial model given 
in Fig. ^ are shown. Here, Qc and pc are, respectively, the 
angular velocity and density at the center of the model. 
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t=2D.n tn 




FIG. 3. The development of the bar mode into the nonhn- 

ear regime is shown using 2-D density contours in the equato- 
rial plane for model L3. The maximum (central) density has 
been normalized to unity at the initial time, and the contour 
levels are at 0.5, 0.05, 0.005, and 0.005. The model rotates in 
the counterclockwise direction. 




Fig. 4 





FIG. 4. Density contours in the equatorial plane for the 
later stages of models Dl (a-c) and D2 (d-f) are shown. The 
contour levels are the same as in Fig. ^ and time is measured 
from the initial moments in the respective simulations. 
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Fig. 5 
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FIG. 5. Same as Fig. ^ for models LI (a-c), L3 (d-f), and 
L2 (g-i). 
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t/t, t/t, 

FIG. 6. Plots of the stability parameter /9 are shown for 
models Dl (dotted line) and D2 (solid line) in frame (a) and 
models LI (dotted line), L3 (solid line) and L2 (dashed-dotted 
line) in frame (b). 
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(a) 




FIG. 7. The position of the system center of mass is shown 
as a function of time, (a) Models Dl (dotted hne) and D2 
(soUd line) (b) Models LI (dotted line) and L3 (solid line). 
Note that the use of 7r-symmctry in model L2 prohibits the 
development of center of mass motion. 
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T / Iq 

FIG. 8. The growth of the amplitudes \ Am \ for m = 1 (dot- 
ted line), m = 2 (thick solid line), m — 3 (dashed lino), and 
m = 4 (thin solid line) is shown. The amplitudes wore calcu- 
lated in the equatorial plane in a ring with radius zu = 0.354 
for runs Dl, D2, LI, and L2, and radius = 0.344 for run 
L3. 
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FIG. 9. The gravitational waveform h+ for an observer lo- 
cated on the axis at 9 — 0, ip — oi & spherical coordinate 
system centered on the source is shown as a function of time 
for all models. The quantity plotted is actually rh+, where 
r is the distance to the source. The quantities h+ and r have 
been normalized to {GM/(?ReY and R^, respectively. 
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t 

FIG. 10. The quantity r/inorm = r {h% + ftx)^^^ is shown 
as a function of time for all models. The normalization is 
the same as in Fig. 9. The absolute scale of the time axis is 
not labeled as the L3 and D2 curves have been horizontally 
shifted in order to align their intial peaks with that of the L2 
curve. 
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Ref. A'' P code 7r-symm tba 



iflnal Atb 



0.33 
0.33 
0.31 
0.30 
0.30 
0.30 
0.30 
1.5 0.327 
1.5 0.304 
1.5 0.327 



1.5 
1.5 
1.8 
1.5 
1.5 
1.5 
1.5 



Eulerian 

SPH 
Eulerian 
Eulerian 
SPH 
SPH 
Eulerian 
Eulerian 
Eulerian 
Eulerian 



yes 
no 
yes 
no 
no 
no 
yes 
no 
no 
no 



2.5 tc 
2.0 tc 
11.3 tc 
10 tc 
8.2 tc 

5.7 tc 

6.8 tc 
6.8 tc 

10.1 tc 

n/a 



9.5 
9.5 tc 

19.4 tc 

15.5 tc 
16 

15.9 
24.3 tc 

12.3 tc 

14.4 tc 
11.2 tc 



> 7.0 tc 

> 7.5 tc 

> 8.1 tc 
~ 5.5 tc 
~ 7.8 tc 
~ 7.9 tc 

> 17.5 tc 

> 5.5 tc 

> 4.3 tc 
n/a 



remarks 
central bar at tflnai 
central bar at tflnai 
central bar at tflnai 
no bar at tflnai 
no bar at tflnai 
bar gone by t ~ 13.6 tc 
central bar at tflnai 
central bar at tflnai 
central bar at tflnai 
central bar at tflnai 



TABLE I. Properties of long-duration simulations of the 
bar mode instability, tbar max is the time at which the bar 
reaches its maximum elongation and tflnai is the end of the 
simulation. The length of time the bar persists is Atbar. In 
this table, time is measured in units of tc, where 1 tc = 1 
central initial rotation period (cirp). 



model code grid size 7r-symmetry piow persistent | Jbar max — Jo|/Jo 
iV^ X iV^ X bar 



Dl 


V 


64 X 64 X 64 


no 


10- 


lU 


no 


2.6 X 10" 


2 


D2 


V 


64 X 64 X 128 


no 


10" 


10 


no 


5.1 X 10" 


2 


LI 


C 


64 X 64 X 128 


no 


10" 


-7 


no 


5.0 X 10" 


3 


L2 


C 


64 X 64 X 64 


yes 


10" 


-7 


yes 


5.0 X 10" 


3 


L3 


C 


128 X 128 X 128 


no 


10" 


-7 


yes 


5.2 X 10" 


3 



TABLE II. Basic properties of the models run on the U 
and C codes are shown. Note that model L2 was actually run 
with TT-symmetry, giving an effective resolution of 128 zones 
in the angular direction. pi o„ is the minimum density in the 
"vacuum" regions; see Sec. IIIB, The last column shows the 



loss of total angular momentum J at the time that (5 reaches 
its first local minimum and the bar reaches it maximum elon- 
gation for each model; cf. Fig. ^. Jo is the initial total angular 
momentum in the model. 



model 


din \A2\/dt 


din \A4,\/dt 


(T2 


0-4 


W2 
















Dl 


0.54 


0.98 


1.5 


3.2 


0.75 


0.80 


D2 


0.55 


1.1 


1.5 


3.0 


0.75 


0.75 


LI 


0.55 


0.94 


2.0 


3.9 


1.0 


0.98 


L2 


0.55 


1.0 


2.0 


4.0 


1.0 


1.0 


L3 


0.55 


1.1 


2.0 


4.0 


1.0 


1.0 



TABLE III. The growth rates, eigenfrequencies, and pat- 
tern speeds for the m = 2 and m = 4 modes are given for all 
runs. Notice that the pattern speeds W2 ~ Wa for all mod- 
els, indicating that the m = 4 mode is a harmonic of the bar 
mode. 
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